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Abstract 



Streamer ionization fronts are pulled fronts propagating into a linearly unstable state; 
the spatial decay of the initial condition of a planar front selects dynamically one specific 
long time attractor out of a continuous family. A stability analysis for perturbations in the 
transverse direction has to take these features into account. In this paper we show how to 
apply the Evans function in a weighted space for this stability analysis. Zeros of the Evans 
function indicate the intersection of the stable and unstable manifolds; they are used to 
determine the eigenvalues. Within this Evans function framework, a numerical dynamical 
systems method for the calculation of the dispersion relation as an eigenvalue problem is 
defined and dispersion curves for different values of the electron diffusion constant and of 
the electric field ahead of the front are derived. Numerical solutions of the initial value 
problem confirm the eigenvalue calculations. The numerical work is complemented with an 
analysis of the Evans function leading to analytical expressions for the dispersion relation 
in the limit of small and large wave numbers. The paper concludes with a fit formula for 
intermediate wave numbers. This empirical fit supports the conjecture that the smallest 
unstable wave length of the Laplacian instability is proportional to the diffusion length 
that characterizes the leading edge of the pulled ionization front. 
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1 Introduction 

1.1 The streamer phenomenon, ionization fronts and Laplacian instability 

A streamer is the first stage of electric breakdown in large volumes, it paves the way of sparks 
and lightning, but also occurs without successive breakdown in phenomena like sprite dis- 
charges above thunderclouds or in corona discharges used in numerous technical applications. 
Recent reviews of relevant phenomena can be found in [20, 45J. Considered as a nonlinear 
phenomenon, the streamer is a finger shaped ionized region that propagates by self generated 
field enhancement at its tip into nonionized media. It has multiple scales as described in [20 ; 
as a consequence one can investigate a hierarchy of models on different levels of refinement 
that are reductions of each other, starting from the reduction from a particle to a continuum 
model |32j to the reduction from a continuum model to a moving boundary model [9] up to 
the formulation of effective models for complete multiple branched streamer trees without 
inner structure that are known as "dielectric breakdown models" \38\ [39j |4"01 [7] . All these re- 
ductions are the subject of current research; the present paper analyzes the stability of fronts 
in the continuum model; the resulting dispersion relation provides a test case for moving 
boundary approximations. 

Specifically, simulations of the simplest continuum model for negative streamers [15 ^ 16 1 147] 
have established the formation of a thin boundary layer around the streamer head. This layer 
is an ionization front that also carries a net negative electric charge. (Positive streamers 
with positive net charge occur as well, but are not the subject of the present study.) The 
configuration of the charge in a thin layer leads to the above mentioned field enhancement at 
the streamer head that creates high ionization rates and electron drift velocities and hence lets 
the streamer rapidly penetrate the non-ionized region. More recent numerical investigations 
show that the boundary layer or front can undergo a Laplacian instability that generates the 
streamer branch [U [42~1 [36l I37j . (We remark that an additional interaction mechanism in 
composite gases like air somewhat modifies this picture |33j while the present analysis applies 
to negative streamers in simple gases like pure nitrogen or argon.) 

1.2 Moving boundary layers and the transversal instability of pulled fronts 

The streamer can be considered as a phenomenon where an ionized phase is separated from 
a non-ionized phase by a moving thin front. This concept [241 0] implies that streamers 
show similar features as moving boundary problems like viscous fingers, solidification fronts 
propagating into undercooled liquids, growth of bacterial colonies or corals in a diffusive field 
of food etc. Quantitative predictions within such models require a proper understanding of 
the front dynamics, in particular, of their stability against perturbations in the transversal 
direction. This stability determines whether perturbations of the front position will grow or 
shrink, and on the long term whether the streamer will branch or not. As a first insight, one 
would therefore like to analyze the stability of planar fronts against transversal perturbations, 
more specifically, the growth or shrinking rate s(k) of a linear perturbation with transversal 
wave length 2ir/k. 

The ionization front in the model for a negative streamer in a pure gas as treated in 
[I3II61I171IM1I251I1I121I361E7], including electron diffusion, creates a so-called pulled 
front that has a number of peculiar mathematical properties: (i) for each velocity v > v * , 
there is a dynamically stable front solution where the stability is conditional on the spatial 
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decay of the perturbation, hence the long time dynamics is selected by the spatial decay 
of the initial front for z — > oo (where z is the spatial variable along the front); (ii) the 
convergence towards this front is algebraically slow in time [21} I22j: (in) this slow dynamics 
is determined in the leading edge of the front that in principle extends up to z —* oo and in the 
dynamically relevant space it will cause Fredholm integrals in the linear stability analysis to 
diverge, therefore curvature corrections cannot be calculated in the established manner [23], 
(iv) the unconventional location of the dynamically relevant region ahead of the front also 
requires particular care in numerical solutions with adaptive grid refinement [37]. For the 
calculation of the dispersion relation, which can be phrased as an eigenvalue problem for s(k) , 
these features pose two challenges: first, the condition on the one-dimensional dynamical 
stability and algebraic convergence properties, which are typical for pulled fronts, will lead 
to an apparently degenerate eigenvalue problem. Second, in a neighborhood of the origin, 
the dispersion curve s(k) is near the continuous spectrum. Hence numerical calculations of 
the eigenvalue problem with finite difference, collocation or spectral methods often lead to 
spurious eigenvalues. A dynamical systems method involving stable and unstable manifolds 
avoids this problem. The stable and unstable manifolds are at least two-dimensional and an 
exterior algebra approach is employed to calculate the manifolds accurately. 

In \17\ 111 l3j. the treatment of pulled fronts and more-dimensional stable/unstable mani- 
folds was circumvented by neglecting the electron diffusion that acts as a singular perturba- 
tion. In this way, the leading edge of the front together with its mathematical challenges is 
removed and the eigenvalue problem can be solved using shooting on the one-dimensional sta- 
ble/unstable manifolds. The resulting problem is characterized by two length scales, namely 
the length scale lixjk of the transversal perturbation, and the longitudinal length scale of 
electric screening through the front that will be denoted by i a . The dispersion relation in 
this case shows a quite unconventional behavior, namely a short wave length instability whose 
consequences are further investigated in [351 IE]- In the present paper, we analyze the dis- 
persion relation including diffusion, mastering the above challenges and deriving quantitative 
results through a combination of analytical and numerical methods. 

1.3 The Evans function and pulled fronts 

The Evans function is an analytic function whose zeros correspond to the eigenvalues of a 
spectral problem, usually a linearization about a coherent structure like a front or solitary 
wave. It was first introduced in [26] and generalized in [I]. In the last decade, the Evans 
function has been applied in the context of many problems and various extensions and gener- 
alizations have been found, see the review papers |30| EH] and references in there. One of the 
first uses of the Evans function in the analysis of a planar front can be found in [46| . which 
analyzes the stability of a planar wave in a reaction diffusion system arising in a combustion 
model. In the current paper we will show how pulled fronts can be analyzed with the Evans 
function by using weighted spaces in its definition. 

To define the Evans function, one writes the eigenvalue problem as a linear, first order, 
dynamical system with respect to the spatial variable z. Along the dispersion curve s(k) , the 
dynamical system has a solution which is bounded for all values of z. This can be phrased 
in a more dynamical way as: the manifold of solutions which are exponentially decaying for 
z — * +oo (stable manifold) and the manifold of solutions which are exponentially decaying 
for z —* — oo (unstable manifold) have a non-trivial intersection along the dispersion curve. 
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The Evans function is a function of the spectral parameters s and k , which vanishes if the 
stable and unstable manifolds have a non-trivial intersection. Hence the Evans function can 
be viewed as a Melnikov function or a Wronskian determinant, see also [29 or references in 
there. 

In case of a pulled front, the definition of the stable manifold, and hence the Evans 
function, is not straightforward. The temporal stability of the asymptotic state of the pulled 
front at +00 is conditional on the spatial decay of the perturbation. So this decay condition 
should be included in the definition of the stable manifold, otherwise the dimension of this 
manifold might be too large. We will show that this condition can be built in the definition 
of the stable manifold by considering the stable manifold in a weighted space. The Evans 
function is defined by using the weighted space for the stable manifold. Hence the dispersion 
curve s(k) can be found as a curve of zeros of this Evans function. 

1.4 Organization of the paper 

In section [21 we recall the model equations and the construction and properties of planar 
fronts. In particular, we summarize the multiplicity, stability, dynamical selection and con- 
vergence rate of these pulled fronts. In section [3j the stability of these fronts is investigated 
as an eigenvalue problem for the dispersion relation s(k) of a linear perturbation with wave 
number k. The dispersion relation depends on the far electric field E^ and the electron 
diffusion D as external parameters. In the stability analysis of the pulled ionization fronts, 
a constraint is imposed on the asymptotic spatial decay rate of the perturbations. This 
constraint corresponds to the decay condition for the one-dimensional stability, but has to 
be chosen quite subtly to avoid problems with the algebraic decay of the front solution. A 
consequence of the decay condition is that the eigenvalue problem (dispersion relation) is 
solved in a weighted space. In this weighted space, the apparent degeneracies have disap- 
peared, the stable and unstable manifolds of the ODE related to the eigenvalue problem are 
well-defined and intersections of those manifolds are determined by using the Evans function. 
In section 13.41 dispersion relations for positive s are derived numerically for a number of 
pairs of external parameters (i?^,!)). The numerical implementation of the Evans function 
uses exterior algebra to reliably solve for the higher dimensional stable and unstable man- 
ifolds. In section [4j the numerical dispersion relation is tested thoroughly and confirmed 
with numerical simulations of the initial value problem for the complete PDE model for 
the particular values (E^^D) = ( — 1,0.1) where D = 0.1 is typically used for nitrogen 
[El[Tl|47l[H[2lll3[ia[Ml[37]and J E; oo = -lisa representative value for the electric field. 
The later sections treat either general (E^^D) analytically or a larger range of {E^^D) 
numerically. 

In section [5j explicit analytical asymptotic relations for the dispersion relation s(k) are 
derived for the limits of small and large wave numbers k . For k = , two explicit eigenfunc- 
tions are known (which are related to the translation and gauge symmetry in the problem). 
These explicit solutions lead to expressions for the solutions on the stable manifold for small 
wave numbers. The interaction between the slow and fast behavior on this manifold leads 
to an asymptotic dispersion relation for small k. For large wave numbers, the eigenvalue 
problem for the dispersion relation is dominated by a constant coefficient eigenvalue problem. 
An eigenvalue exists only if this constant coefficient system has no spectral gap. Using expo- 
nential dichotomies and the roughness theorem, the asymptotics of the dispersion relation is 
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derived by a contradiction argument. In section [61 these asymptotic limits are tested on the 
numerical data derived in section [3j It is found that the asymptotic limit for small k fits the 
data very well, while the asymptotic limit for large k is not yet applicable in the range where 
s(k) is positive. After a discussion of relevant physical scales, we suggest a fit formula joining 
the analytical small k asymptotic limit with our physically motivated guess. This formula 
fits the numerical data well for practical purposes and strongly supports the conjecture that 
the smallest unstable wave length is proportional to the diffusion length that determines the 
leading edge of the pulled front. 

2 The streamer model and its ionization fronts 

In this section we describe the streamer model and summarize the features of planar ionization 
fronts as solutions of the purely one-dimensional model as a preparation for the stability 
analysis in the dimensions transversal to the front. In particular, we recall the multiplicity of 
the front solutions that penetrate a linearly dynamically unstable state, and the dynamical 
selection of the pulled front. 

2.1 Model equations 

We investigate negative fronts within the minimal streamer model, i.e., within a "fluid ap- 
proximation" with local field-dependent impact ionization reaction in a non-attaching gas like 
argon or nitrogen [241 l25l [T71 SJ H2] . The equations for this model in dimensionless quantities 
are 

d t u - DV 2 a - V-(ctE) = <r/(|E|), (2.1) 

d t p = <7/(|E|) , (2.2) 
V-E = p-a, E = -V<p, (2.3) 

where a is the electron and p the ion density, E is the electric field and 4> is the electrostatic 
potential. For physical parameters and dimensional analysis, we refer to discussions in \24\ 
[25j [T71 IH E2] • Tn e electron current is approximated by diffusion and advection — DVa — crE. 
The ion current is neglected, because the front dynamics takes place on the fast time scale 
of the electrons and the ion mobility is much smaller. Electron-ion pairs are assumed to 
be generated with rate a/(|E|) = o"|E|a(|E|) where cr|E| is the absolute value of electron 
drift current and ck(|E| ) the effective impact ionization cross section within a field E. Hence 
/(|E|) is 

/(|E|) = |E|a(|E|) . (2.4) 

For numerical calculations, we use the Townsend approximation a(|E|) = e -1 /' E ' [241 1251 \T7\ 
HK32]. For analytical calculations, an arbitrary function a(|E|) > can be chosen where we 
assume that a(0) = and therefore /(0) = = /'(0). We will furthermore assume that 
q(|E|) is monotonically increasing in |E|, this is a sufficient criterion for the front to be a 
pulled one [22]. The electric field can be calculated in electrostatic approximation E = —Vcj). 

Mathematically, the model (|2.ip - (|2.3p describes the dynamics of the three scalar fields a , 
p and (p. It is a set of reaction-advection-diffusion equations for the charged species a and 
p coupled nonlinearly to the Poisson equation of electrostatics. 
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2.2 Two types of stationary states 

It follows immediately from (j2.1[ )- (12.3[) that there can be two types of stationary states of the 
system, one characterized by a = and the other by E = (as /(|E|) = implies |E| = 0.). 

The stationary state with a = is the non-ionized state. As the dynamics is only carried 
by the electrons a, there is no temporal evolution for a = even if the ion density p has 
an arbitrary spatial distribution. The electric field E = — V</> then is determined by the 
solution of the Poisson equation — V 2 (/> = p and by the boundary conditions on 0. In certain 
ionization fronts in semiconductor devices [43 J , it is essential that the equivalent of p does not 
vanish in the non-ionized region. In the gas discharges considered here, on the other hand, it 
is reasonable to assume that the non-ionized initial state with a = also has a vanishing ion 
density p = 0, and therefore no space charges. 

The stationary state with vanishing electric field E = describes the ionized, electrically 
screened charge neutral plasma region behind an ionization front, the interior of the streamer. 
From E = the identity V • E = follows immediately, and therefore electron and ion 
densities have to be equal a = p. In the absence of a field, the electrons diffuse dto = DV 2 a 
while the ions stay put dtp = 0. Therefore, these densities only can stay equal if V 2 p = 0. 
Simulations [IH HH [47l Ell EH HI S21 [36l E7] show that this occurs typically only if p is 
homogeneous (though counter examples can be constructed). 

2.3 Planar ionization front solutions 

An ionization front separates such different outer regions: an electron-free and non-conducting 
state with an arbitrary electric field ahead of the front from an ionized and electrically 
screened state with arbitrary, but equal density a~ = p~ of electrons and ions. In particular, 
we are interested in almost planar fronts propagating into a particle-free region p = a = 
(where therefore V 2 = 0), and we study negative fronts, i.e., fronts with an electron surplus 
that propagate into the electron drift direction towards an asymptotic electric field Eoo < . 
For a planar front, it follows from V 2 ^ = — V • E = that the electric field ahead of the front 
is homogeneous. 

We assume that the front propagates into the positive z direction; the electric field ahead 
of a negative front then is E — > EaJz , < , for z —* oo . (Here z is the unit vector in the 
z -direction.) It is convenient to introduce the coordinate system (x,y,£ = z — vt) moving 
with the front velocity v = vz. A planar, uniformly translating front is a stationary solution 
in this co-moving frame, hence it depends only on the co-moving coordinate £, and will be 
denoted by a lower index 0. A front satisfies 

Dd^ao + (v - <%0o) <%o- + o"o(po - 0o) + co/o = 0, 

vd^po + oo/o = 0, (2.5) 

<9|0o + Po ~ co = 0, 

where /o = J(\Eq\). This system can be reduced to 3 first order ordinary differential equa- 
tions. First, due to electric gauge invariance, the system does not depend on 4>q explicitly, but 
only on Eq = —d^(f)Q. Using the variable Eq instead of (f>o reduces the number of derivatives 
by one. Second, electric charge conservation dtq + V • j = can be rewritten in co-moving 
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coordinates for a uniformly translating front as —vd^qo + d^jo = . Therefore it can be inte- 
grated once — vqo + jo = c, <9gc = 0. In the present problem, the space charge is qo = po~ &o 
and the electric current is jo = —Dd^dQ — (JqEq . Furthermore, as there is a region with 
vanishing densities do = = po ahead of the front, the integration constant c vanishes in 
this region, and therefore everywhere. Thus the planar front equations (|2.5p can be written 
as 

Dd^a = v(po- cr ) - E a , 

vd^po = -a f(\E \), (2.6) 
d{E = ,oo-0"o, 

where d^cfro = —Eq decouples from the other equations. The planar front equations imply 
that E (Q < for all £ when E^ < 
The fronts connect the states 



and po - [a- , (2.7) 
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and the electrostatic potential (fro connects <fr~ (for £ — ► — oo) with —E^ + cfr + (for £ — > 
+oo). The ionization density a~ behind the front and the electrostatic potential difference 
(fr + — <fr~ have to be determined for arbitrarily chosen electric field E^ ahead of the front 
and for arbitrary, but sufficiently large, front velocity v . (We remark that only the potential 
difference <fr + — (fr~ matters due to the gauge invariance of the electrostatic potential as one 
easily verifies on the equations.) The fronts can be constructed as heteroclinic orbits in a 
three-dimensional space as demonstrated in [25J. 

The diffusion constant D is obviously a singular perturbation. For D = 0, the front 
equations can be solved analytically |25|, [3], i.e., one can find explicit expressions for the 
particle densities a"o[-Eb]> Po[^o] and for the front coordinate £[Eo] as a function of the 
electric field Eq . For the negative fronts treated here, the front is continuous as function of 
D and the limit D — » exists and equals the value of the front at D = , while for positive 
fronts (Eqo > 0), it is singular 



2.4 Multiplicity of front solutions, pulled fronts and dynamical selection 

The non-ionized state (a,p,E) = (0, 0, £00) with a nonvanishing electric field E^ is linearly 
unstable under the temporal dynamics of the PDE (|2.ip - (|2.3p . In fact, this spatial region 
ahead of the front dominates the dynamics, cf. the discussion in [25, 22J. Therefore, for 
fixed £"00 , there is a continuous family of uniformly translating solutions, parametrized by 
the velocity v > v* [Ml E3 EH [22] , where 



v*^) = |£oc| + 2V-D /(|#oo|)- (2.8) 

The dynamics of uniformly translating fronts with velocity v > v* are dominated by a flat 
spatial profile in the leading edge of the front 



*,(e)^°e-* withA<A* = A /^^, (2.9) 
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where velocity v and decay rate A are related through 

v(E 00 ,\) = \E 00 \+D\ + ^^, (2.10) 

A 

and therefore v(E 00 ,X) > v*{E 00 ) = v(E 00 ,A*) for A / A* . The spatial profile (|2l)|) with 
A < A* cannot build up dynamically from some initial condition with larger A; and it will 
destabilize if perturbed with an initial condition with smaller A , therefore such flat and fast 
fronts can be approached dynamically only by initial conditions with exactly the same profile 
(|2.9p in the leading edge. For a thorough discussion of this dynamics, we refer to [22J. 

In practice, the continuum approximation for the electron density breaks down for very 
small densities in the leading edge and the initial electron distribution satisfies a decay con- 
dition of the form 

lim cr(x,y,£,t = 0) e A? = for all A < A*, (2.11) 

if the penetrated state is really non-ionized. Therefore the velocity v* is called the "selected" 
one, because it is the generic attractor for most physical initial conditions. Mathematically 
speaking, the profile with velocity v* (the selected front) is the only profile that can build up 
dynamically from steeper initial conditions. 

Therefore the condition (|2.1ip on the spatial decay of the initial electron distribution 
excludes all front solutions with velocity v > v* as long time attractors of the dynamics. If 
the criterion (|2.1ip is satisfied, then the selected front with speed v* is dynamically stable 
and is approached with the universal algebraic convergence rate in time |21| [22] 

= „* _ _p_ + O [ JLV (2.12) 
v ; 2A*t \t 3 / 2 J 

However, without the spatial decay condition on the initial condition, the selected front is 
formally not stable (although this is physically irrelevant). This will lead to specific problems 
and solutions in the transverse stability analysis presented in the next section. 
The spatial profile of the electron distribution in the selected front is 

<7 V * (£) 5 ~°° « + b) e" A * ? , a > 0. (2.13) 

To summarize, if the analysis is restricted to initial conditions with a sufficiently rapid spatial 
decay in the electron densities (|2,lip . then the fronts have only one free external parameter, 
namely the field Eoo ; it determines the asymptotic front velocity (|2.8j) and profile fl2. 13f> after 
sufficiently long times. Furthermore, the equivariance in the system gives that the position 
of the front and its electrostatic potential are free internal parameters. 

2.5 Full spatial profiles of the selected pulled planar front 

The spatial decay behind the front will be important in the analysis as well, therefore we 
recall the basic behavior. For £ — * — oo, the electron density approaches 

(0 ^=°° °~ + c e x ^, c > 0, (2.14) 
and the electric field decays with the same exponent E(£) = — (c/A~) e x For D = 0, 

a'(E oo ,D = 0)= a(x)dx (2.15) 

Jo 
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was derived in [25] . For D > , a decreases by a correction of order of D , more precisely, 
a-(E oo ,D) = a-(E oo ,0) + O{D), a~(E oa ,D > 0) < a^(E oo ,0) (2.16) 
was proved in the appendix of [32]. The eigenvalue A - is given by 

vV 2 + ADa~ - v* 



A" 



2D 



(2.17) 



where both v* and a depend on E^ and D . For small D , A can be approximated as 

^ E °°\ a(x) dx 



A" 



a 



+ 0(D) 



< v ./O I -^o 



+ 0(VD) 



(2.18) 



As an illustration, the spatial profiles of electron and ion density and the electric field of 
the selected front solution for a range of fields E^ and diffusion constants D are plotted in 
Figure [TJ 
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Figure 1: The pulled planar front solutions on the left for varying E^ = — 1 , —5 and —10 and 
fixed D = 0.1, and on the right for fixed E^ = — 1 and varying D = 0.1, 0.01 and 0. The 
upper panels show scaled electron and ion densities ao(^)/a~(E (X1 ,D) and po(^)/a~ (Eqo, D) 
and the lower panels the corresponding scaled electric fields .E7o(£)/l-^oo| as a function of 
the spatial coordinate £. The fronts are displayed in a staggered way. The normalization 
factors a~ (Eqc, D) in the upper panels are cr~(— 1, 0) = 0.149, a~(—l, 0.01) = 0.148, 
o-(-l, 0.1) = 0.144, o-(-5, 0.1) = 2.832, cr~(-10, 0.1) = 7.169. 
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3 Numerical calculation of the dispersion relation 

First we will introduce the transversal perturbation setting and discuss an apparent degener- 
acy of the dispersion relation. However, it turns out that the constraint on the spatial decay of 
the electron density "selects" a single dispersion relation for every far field . This relation 
then is calculated numerically based on dynamical systems techniques involving intersections 
of stable and unstable manifolds. Results for different fields E^ and diffusion constants D 
are presented. 



3.1 Linear transversal perturbations of planar fronts 

Suppose that there is a linear transversal perturbation of the uniformly translating front 

<r{x,y,S,t) = a (0+6a 1 {x,y,£,t) + O(5 2 ), £ = z-v% (3.1) 

and similarly for p and <fi. The linearized equation for a\ , ~p± and 4>\ follows from Eqs. (|2. lj) - 
(|2.3p . By decomposing the perturbations into Fourier modes in the transversal directions x 
and y , by using isotropy in the transversal (x, y) plane and by using a Laplace transformation 
in t, the ansatz 



ikx+st 



{<Tk>Pk,<j>k){Q 



(3.2) 



can be used for each Fourier component. The resulting equation can be written as a linear 
first order system of ODEs, using the extra variables = d^ak and E^ = —d^cftk . Introduce 
w = (rfc, cifc, pk, Ek, 4>k) and the linear system is given by 



<%w = M(£;#oo,M) w, 



(3.3) 



with M 



E +v* 
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d v f{rj)\ v=m are used. 



In the matrix M, the abbreviated notations /o = f(\Eo\) and /q 
For the terms with /q , we have used that Eq < , hence = — 1 . 

As the matrix M depends on k 2 , but not on k itself, the matrix is invariant under the 
transformation k — * —k. Thus if s(k) = s* , then also s(—k) = s* and vice versa. Therefore 
it is sufficient to determine the dispersion relation for k > and this will imply the relation 
for k < and from now on, we will use the convention that k > . Note that the invariance 
implies only that the dispersion relation will be a function of \k\. As will be shown later, the 
dispersion relation is not an analytic function of k near k = and its expansion near k = 
is linear in \k\. 

For future use, we remark that the linearization matrix M does not involve any £- 
dependent terms in the fourth and fifth row and implies that Ek and (j>k are related by 



EL 



Thus the E^ -component of any solution of the linearized system (|3.3p can be 
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expressed as an integral 



E k (0 = Cl e^ + c 2 e-^ + 



-fcf 



1 /* 



[pk(v) -°k(v)] dv, 



(3.4) 



where the constants c\ and c% are determined by the value of E^ and at £ = £q • 



Stable and unstable manifolds and degeneracy of the dispersion rela- 
tion 



3.2 



The linearized problem (|3.3|) is a spectral problem with the spectral parameters s and k . If 
the asymptotic matrices M ± (£/ 00 , fc, s) = lim^-too M(£; E^,, fc, s) exist and are hyperbolic 



(i.e., no eigenvalues on the imaginary axis), then the system (|3.3p has a bounded solution 
if and only if the unstable manifold from £ = — oo and the stable manifold from £ = oo 
have a non-trivial intersection. So we will focus in this section on determining the stable and 
unstable manifolds. 

The behavior of the unstable manifold at the back of the front is given by the asymptotic 
matrix 



M-(E 00 ,k,s) = lim M(£; #oo,M) 



For s > and k ^ 0, this matrix has two negative and three positive eigenvalues: 



I D 


cr-+s+Dk 2 
D 


(7 

D 





o \ 


1 




















s 

V* 











-1 


1 





-k 2 


V o 








-1 


o ) 



V 



yV 2 + 4D(<T- + S + Dfc 2 ) 

ZD ' 







(3.5) 

fe) is identical 



Thus the unstable manifold is three dimensional. We remark that /uT(s 
to the spatial decay rate A - (|2.17p behind the unperturbed front. 

Finding the stable manifold ahead of the front is less straightforward. Normally the stable 
manifold ahead of the front would be characterized by the matrix liiri£_> +00 ~M.(£ t ; E^, k, s) . 
For s > and s + Dk 2 < f{E 00 ) this matrix exists and has two positive and three negative 
eigenvalues: 

I s + Dk? _ -y/f(E^j±Vs + Dk 2 



v* 



IB 



(3.6) 



Thus the stable manifold is three dimensional and a dimension count gives that the inter- 
section of stable and unstable manifold is generically one dimensional. So for small values 
of s and k, a continuous family of eigenvalues seems to exist. This feature is related to the 
instability of the asymptotic state at +oo, to the continuous family of uniformly translat- 
ing solutions for all v > v*(E 00 ), and to the instability of fronts against perturbations with 
smaller spatial decay rate A, as discussed in the previous section. The continuous family of 
eigenvalues s for fixed wave number k is eliminated by applying the analysis only to fronts 
with a sufficiently rapid spatial decay (12. lip . This condition will be imposed in the definition 
of the stable manifold; it will make the spectrum discrete. 
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Define the scaled vector 



w = D w, D = diag(e (A *-^ , e ( A *-^, 1, 1, 1) 



(3.7) 



where (3 6 (0, A*) will be fixed later and depend on k and A*. The freedom in the choice 
of (3 is reminiscent of the fact that the decay condition holds for any A < A* , but not for 
A = A* . The system for w is 

w s = M(£;£oo,Mw, (3.8) 

with 



M = D M D 1 + (dtp) D 



D 



+ A* - (3 
1 


2a - Po -fo+s+Dk 2 


_£0 JA*-/3)t 
D e 






\ 



D 

A* -(3 


D e 







-fo p -(A*-/3)£ 
v* c 

_ e -(A*-0)£ 


s 

u* 
1 







-k 2 











-1 


o / 



V 

Note that if = , then the asymptotic matrix ahead of the front (at £ = +oo ) does not 
exist because e A *£<7o(£) grows linearly in £ according to (|2. 13|) . To get an asymptotic matrix 
ahead of the front, it is necessary that < (3 < A* . In this case, the asymptotic matrix is 



lim M(£;Eoo,fc,s) 

§— >oo 



/- 


-A* - [3 


-foc+S + Dk 2 








\ 


D 




1 


A* -(3 
















foe 


s 

V* 






















-k 2 


V 











-1 


/ 



where = /(|-Eoo|). The matrix M + has the eigenvalues 



v 



and 



s + Dk 2 



D 



(3. 



Hence for s > and < (3 < min(A*, k\J\ + s/(Dk 2 )) , there are two negative and three 
positive eigenvalues. Thus the stable manifold of (|3.8p is two dimensional. 

For the original unsealed system (|3,3p this means that only the two-dimensional sub- 
manifold given by D _1 acting on the stable manifold of (|3.8|) is relevant for the transverse 
instability. This submanifold will be called the stable manifold of (|3.3p from now on. With 
this convention, the dispersion relation is a well-defined curve s(k) and the curve is such that 
at s = s(k), the linearized system (13. 3p has a bounded solution which satisfies the spatial 
decay condition (|2.1ip . Note that for both asymptotic matrices M + and M_ , the eigenvalues 
ztk become a degenerate eigenvalue at k = 0. This leads to square root singularities and 
it can be expected that the dispersion relation s(k) will be a function of y/W = \k\ for k is 
small. This will be confirmed in section [SJ 
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3.3 The Evans function for the transverse stability problem 

The occurrence of an intersection of the stable and unstable manifolds will be measured 
with the Evans function. Our numerical method to determine the dispersion curve as an 
eigenvalue problem is based on a definition of the Evans function in an exterior algebra 
framework and uses similar ideas as in [21 \TT\ [T2l [8j [TP] [T4"] . The approach of following the 
stable/unstable manifolds at £ = ±oo with a standard shooting method and checking their 
intersection using the Evans function, works only if these manifolds are one-dimensional or 
have co-dimension one; in the present model, this is the case in the singular limit D = and 
a shooting method was used in [3]. Otherwise, any integration scheme will inevitably just 
be attracted by the eigendirection corresponding to the most unstable (stable) eigenvalue. 
Exterior algebra can be used to avoid this problem for higher dimensional manifolds and 
to preserve the analytic properties of the Evans function. Recently, a different method to 
calculate the Evans function for higher dimensional manifolds has been proposed in [27J. 
This method uses a polar coordinate approach and looks like a more suitable method for very 
high dimensional problems. 

To calculate the evolution of the two dimensional stable and three dimensional unstable 
manifold in a reliable numerical way, we will use the exterior algebra spaces A 2 (C 5 ) and 
/\ 3 (C 5 ), respectively. The advantage of these spaces is that in /\ (C ra ), an I -dimensional 
linear subspace of C n can be described as a one-dimensional object, being the I -wedge product 
of a basis of this space. Also, the differential equation on M 5 (or C 5 ) induces a differential 
equation on the spaces f\ (C 5 ) : 

Wt = M®fcE 00 ,k,8)W, WeA'(C 5 ). (3.10) 

Here the linear operator (matrix) is defined on a decomposable /-form wi A ... A wj, 

Wj G C 5 , as 

M W (wi A ... A wj) := (Mwi) A...Awj + ...+Wi A...A (Mw ( ) (3.11) 

and it extends by linearity to the non-decomposable elements in /\'(C 5 ). General aspects 
of the numerical implementation of this theory can be found in [2\. The general form of the 
matrices M( 2 ) and M( 3 ) can be found in the appendix. 

To determine the three-dimensional unstable manifold for £ G (— oo, 0] , we will use (|3.10p 
with 1 = 3. Since the induced matrix M^ 3 ) (£; E^, k, s) inherits the differentiability and 
analyticity of M(£; E^, k, s) , the following limiting matrix exists: 

M®(E 0O ,k,s) = lim M^; 

oo 

The set of eigenvalues of the matrix Mj_ ; (£' 00 , k, s) consists of all possible sums of three 
eigenvalues of 'M.±(E O0 , k, s) (see Marcus [34]). Therefore, for s > and k ^ 0, there is an 

(3) 

eigenvalue za_ of M_ , which is the sum of the 3 positive eigenvalues of ML , i.e., 

s v* Jv* 2 + 4D(cr- + s + Dk 2 ) 

V- = k -\ — + — — 

v* 2D 2D 

(note that the subscript " — " in v- refers to exponentially decaying behavior at — oo, not to 
the sign of z/_, which is obviously positive). The eigenvalue z/_ is simple and has real part 
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strictly greater than any other eigenvalue of M_ (as ML is hyperbolic). We denote the 

— (3) — — 

eigenvector associated with as W~ , i.e., Ml'W^ = i/_W" . This vector can always be 
constructed in an analytic way (see [31j pp. 99-101], [101 HI HH]). In this case it is easy to 
determine an explicit analytical expression for the eigenvector as M_ is quite sparse. The 
unstable manifold corresponds to the solution W - (£) of the linearized system (|3.10p (with 
1 = 3) which satisfies lim e~ u ~^~W~ (£) = W~ . 

>— oo 

The stable manifold can be determined in a similar way. As indicated in the previous 
section, the scaled system (|3.8p will be used to determine the stable manifold. For the stable 
manifold with £ G [0, oo) , we will use (|3.10j) with I = 2 and the scaled matrix M . As before, 
the asymptotic matrix 

M^ 2) (#oo,M = lim M(%;£oo,M. 

(2) 

exists. Now the eigenvalues of M\ (-#005 s) consists of all possible sums of two eigenvalues 

(2) 

of M±(E 00 ,k,s). Therefore, for s > 0, k^O, M K +> has an eigenvalue v + , which is the sum 
of the 2 negative eigenvalues of M+ , i.e., 

/ Is + Dk 2 \ 
"+ = -(V-D- + *-' 3 j 

As before, this eigenvalue is simple and has real part strictly less than any other eigenvalue of 

(2) i (2} 

. The eigenvector associated with u + will be denoted by Wj , i.e., = ^ + W+ 

The stable manifold of the scaled system (|3.8[) corresponds to the solution W + (£) of the 
linearized system (f3T0|) (with I = 2 and M = M) which satisfies lim e^+^W+(£) = W+ . 

To get the stable manifold of the original unsealed system, the inverse scalings matrix D -1 (£) 
has to be used. For arbitrary £ > , the transformation in the wedge space A 2 (C 5 ) is quite 
complicated, but we will only need the original stable manifold at £ = 0. And at £ = 0, the 
scalings matrix is the identity matrix. Hence at £ = 0, the scaled stable manifold and the 
original stable manifold are the same and W+ (0) describes the stable manifold of (|3.3p at 

e = o. 

With the stable and unstable manifold as found above, the Evans function can be defined 

as 

A(E 00 , k, s) = W~(0; E^, k, s) A W + (0; -Eoo, s), s>0,k^0. (3.12) 

Thus the Evans function A is more or less the determinant of the matrix formed by a basis 
of the unstable manifold at £ = and a basis of the stable manifold at £ = . If this function 
is zero, then the bases are linearly dependent, hence the two manifolds have a non-trivial 
intersection. 

We have focused on the case s > 0. For —Dk 2 < s < 0, the system is still hyperbolic, 
but with a two dimensional unstable manifold and a three dimensional stable manifold. The 
method above can be easily adapted to calculate the dispersion curve in this region too. 

3.4 Numerical results on the dispersion relation with the Evans function 

To calculate the Evans function numerically, first the front solution has to be determined 
numerically as it appears explicitly in the linearization matrix M(£; E^, k, s) . The front is an 
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invariant manifold connecting two fixed points of the ODE (|2.6h . so it can be easily determined 
by invariant manifold techniques or shooting, using the package DSTool [6j. Shooting works 
in this case as the front connects a one-dimensional unstable manifold to a three-dimensional 
center-stable manifold in the ODE f)2.6|) . 

After determining the fronts, the stable and unstable manifolds can be calculated by 
numerical integration, see e.g. [2j 1101 [12], In the numerical calculation of the stable manifold, 
we will use = \ min(A*, k) . For the stable manifold, the linearized equation on /\ 2 (C 5 ) 



M( 2 )(£;£oo,M-M£oo,Mll W+ W + (£)L, = W+(E 00 ,k,s) 



\£=L 

is integrated from x = to £ = 0, using the second order Gauss-Legendre Runge-Kutta 
(GLRK) method, i.e. the implicit midpoint rule. Here the scaling W + (£) = e~ v+ ^ W + (£) en- 
sures that any numerical errors due to the exponential growth are removed and W + (£) | = 
W + (^)|^_ Q is bounded. The eigenvector W+(£' 00 , k, s) can be determined explicitly as wedge 
product of the relevant eigenvectors of M + (i? 00 , s, k) thanks to the sparse nature of this ma- 
trix. 

For the unstable manifold, the linearized equation on /\ 3 (C 5 ) 



w- 



M^;£oo,M)-M£oo,*,fc)ll W", W-(0L =r =W~(E 00 ,s,k) 



\£=L B 

is integrated from x = —-Loo to £ = , also using the implicit midpoint rule and introducing 
the rescaling W~(£) = e' u ~^ W~(f) to remove potential exponential growth. Again, the 
eigenvector W~(^ 00 ,fc, s) can be determined explicitly as wedge product of the relevant 
eigenvectors of M.~(E oc , s, k) . 

At £ = 0, the computed Evans function is (see (|3.12p ) 

A(S 00 ,fc,s) = W~(0) AW + (0) = W~(0) AW+(0). (3.13) 

For s = = k , the center-stable and the center-unstable manifold have a two-dimensional 
intersection, due to the translation and gauge invariance, see section 15.11 for details. In order 
to determine the dispersion curve, we start near k = and s = and then slowly increase k 
and determine for which s(k) the Evans function A(^ QO , k, s(k)) vanishes. 

The numerical errors in the calculation of the Evans function are mainly influenced by 
the step size used in the numerical integration with the GLRK method and errors in the 
numerically determined front. The numerical integration uses the step size Sx = 0.01. We 
have performed various checks with a decreased step size and this shows that the error in the 
value of s for fixed k is largest (order 10 -4 ) if k is small and decreases for larger k (order 
10~ 6 ). The accuracy of the front has been checked and is such that the error in the front 
gives a negligible error (compared to the error due to the error in the step size) in the value 
of s(k) . It turns out that the scheme is not very sensitive to errors in the front (at least for 
the .Eoo and D values considered). 

In the following sections, we will present data for the dispersion curve for varying electric 
field and diffusion coefficient D . A more detailed discussion of the data, relation with 
analytical asymptotics and some empirical fitting can be found in section El 

3.4.1 Varying the electric field ahead of the front 

First we consider how the dispersion curve depends on the electric field Eoo ahead of the front, 



while the diffusion coefficient is fixed to D = 0.1 . In Figure 2(a) , the dispersion curve is shown 
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for = — 1, Ecq = —5 and E^ = — 10. The figure shows that the shape of the dispersion 
curve stays similar, but the scales of s and k increase when E^ increases. The dispersion 
curves can be characterized by the maximal growth rate s max and the corresponding wave 
number fc max where s(/c max ) = s m ax as well as by the wave number ko > with s(ko) = 
that limits the band < k < ko of wave numbers with positive growth rates. 




Figure 2: Dispersion curves s(k): (a) for varying E^ and fixed D = 0.1, and (b) for fixed 
Eoq = — 1 and varying D . The pairs (E^, D) shown are the same as in Fig. [TJ The data for 
the singular limit D = is taken from [3]. 



3.4.2 Varying the diffusion coefficient 

Next we consider the effect of varying the diffusion coefficient D, while keeping the electric 
field ahead of the front fixed at E^ = — 1 . In |3| it is shown that if diffusion is ignored 
(D = 0), the dispersion curve stays positive and is monotonically increasing to the saturation 
value s(k) = f(\E OD \)/2 for k — » oo. Our numerics shows that if diffusion is present, this 
is not the case anymore. This is not surprising as diffusion is a singular perturbation. In 



Figure 2(b) the dispersion curve is shown for D = 0.1, D = 0.01 and D = 0; the data for 
D = is taken from [3]. It shows that the growth rate s(k) has a maximum s max if diffusion 
is present and becomes negative for k larger than some ko . Furthermore for decreasing 
diffusion D , the maximal growth rate moves upward towards the saturation value /(l^ooD/2 
for D = 0. This suggests that some features of the dispersion curve behave regularly in D, 
in spite of the fact that D is a singular perturbation. For example, for a finite wave number 
interval, the limit of the dispersion curves for D — > exists and is the curve for D = 0. 
However, the asymptotic profile for large values of the wave number is obviously singular 
in D. This duality can also be found in the front itself: the velocity and the profile of the 
ionization density and the electric field of the uniformly translating negative front depend 
regularly on D = 0, while the profile of the ionization density is singular, as discussed in 
section 12.41 and shown in Fig. [TJ 
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4 Numerical simulation of the perturbed initial value problem 

In the previous section, we have determined the dispersion relation s(k) for transversal pertur- 
bations of ionization fronts as a temporal eigenvalue problem of the PDE system linearized 
about the uniformly translating planar front. Since we are dealing with pulled fronts (cf. 
sections [T1 and [2~4"|) . the problem is unconventional: both the velocity v* of the uniformly 
translating planar front and the dispersion relation s(k) of its transversal perturbations are 
unique only if the spatial decay constraint (|2.1ip is imposed. Furthermore a longitudinally 
perturbed planar front approaches its asymptotic profile and velocity algebraically slowly in 
time (|2.12p . Therefore it is worthwhile to test the predicted dispersion relation on direct 
numerical simulations of the corresponding initial value problem. 

In this section, we will therefore simulate the temporal evolution of a perturbed planar 
front by numerically solving the full nonlinear PDEs (|2.1|) - (|2.3p . and we will determine the dis- 
persion curve from a number of simulations with perturbations with different wave vectors k . 
This is done for far field = — 1 and diffusion constant D = 0.1. 

To determine the instability curve with a simulation of the full PDE, we parametrize the 
evolution of a perturbed planar front with wave number k as 

XJ(x,z,t)^Uo(ti)+5U 1 (ti,t)e lkx+st , £ = z-v*t, V=(a,p,<f>). (4.1) 

If 5 e st is small enough, the solution is in the linear regime, and s(k) can be determined 
from the evolution of the perturbation after Ui(£, t) has relaxed to some time independent 
function. Therefore in the numerical simulations, we choose 5 for each wave number k in 
such a manner that t is large enough to extract meaningful growth rates and that 5 e st is 
small enough that the dynamics at the final time is still well approximated by the linearization 
about the planar front. 

Furthermore, an appropriate choice of the initial condition reduces the initial transient 
time during which Ui(£, t) in the co- moving frame still explicitly depends on time t. Ideally, 
such an initial condition is of the form XJ(x, z, 0) = Uo(£) + 5Ui(£) coskx etc., where Ui 
is a solution of the linearized system (|3.3[) . To find an approximation for Ui(£), we use 
that the instability acts on the position of the front, i.e., we write the perturbed front as 
U {£ + 5e ikx+st ) » U (£) + 5 e ikx+st d c U (C) . Therefore we choose Ui(f) = 8fcU (f) and the 
initial condition as 

U(x,z,0) = U (z) + 5d z U (z) coskx. (4.2) 

As <9^Uo(^) is a solution of the linearized system for k = = s , this choice will be very 
efficient for small values of k and require longer transient times for larger k. 

To solve the full 2D PDE, the algorithm as described in [H [32] is used, while adaptive 
grid refinement as introduced in [37] was not required. For fixed k, the PDE with initial 
condition (|4.2j) is solved on the spatial rectangle (x, z) G [0, L x ] x [0, L z \ . The length of the 
domain in the transversal x -direction, L x , is such that exactly 5 wave lengths fit into the 
domain, i.e., L x = and periodic boundary conditions are imposed in this direction by 
identifying x = with x = L x . On the boundaries in the longitudinal z -direction, Neumann 
conditions for the electron density are imposed. The potential is constant far behind the 
front and the electric field is constant far ahead of the front; therefore for the potential 4>, 
the Dirichlet condition <p = is imposed at z = 0, and the Neumann condition d z (j) = — -Eoo 
at z = L z accounts for the far field ahead of the front. 
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The amplitude of the perturbation is conveniently traced by the maximum of the electron 
density 

Pmax(x,t) = max cr(x,z,t) (4.3) 

evaluated across the front. The reason is as follows. First, Figure [1] shows the spatial profiles 
of planar fronts for different electric fields and illustrates that for fixed D , the maximum 
of the electron density <r max as well as the asymptotic density o~ behind the front strongly 
depend on the field E + immediately ahead of the front, where for a planar front the close 
and the far field are identical: E + = E^ . Second, the modulation of the front position 
leads to a modulation of the electric field E + immediately before the front (cf. discussion in 
section HT2j) : therefore <7 max as a function of E + is modulated as well. 



- + H — I 



+ + + 



+ + 



(a) The maximal value of the electron density 
ffmax(i,t) for t — 50 as a function of the transversal 
coordinate x . The perturbation has wave number 
k — 0.45, the transversal length L x = Wn/k leaves 
space for 5 wave lengths that are clearly visible. 



(b) The logarithm of the amplitude of the front mod- 
ulation log A as a function of time t for the same k . 



Figure 3: Examples of data of the initial value simulation from which the growth rate s(k) 
shown in Fig. 0] are determined. 

An example of <r max (a;,t) as a function of the transversal coordinate x for a fixed time t 



is plotted in Fig. 3(a) . The amplitude of the wave modulation is determined by the Fourier 
integral 

k 



A(t, k) 



10- 

k 



(x, t) cos kx dx. 



In Figure [3(b)| we plot log A against time t for k = 0.45. Note that k = 0.45 is close to 
fop = 0.482 (see Figure 2(a) and Table[T|) where the growth rate vanishes, s(ko) = 0, therefore 



the growth rate in the present example is small and particularly sensitive to numerical errors. 

Figure [3 (b) | shows an initial temporal transient before steady exponential growth is reached 
(where exponential growth manifests itself as a straight line in the logarithmic plot). This is 
typically observed for the larger k -values (k > 0.1); as said before, this is related to the fact 
that the function Ui(z) in the initial condition (14. 2j) is not optimal. For k < 0.1, there are 
less transients as Ui(z) « d z XJo(z) for small values of k. 
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To determine the growth rate s(k) , a least squares algorithm is used to fit the best line 
through the data points (t, log ^4) , and the initial transient time is ignored for larger values 
of k. For each value of k, the growth rate is determined with various choices of 5. The 
resulting growth rate s(k) is indicated in Figure [4] with crosses X and the error bars are 
related to the various choices of 5. 
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Figure 4: The dispersion curve s(k) for = —1 and D = 0.1. The crosses x with error 
bars indicate results of simulations of the full initial value problem as discussed in section U] 
and demonstrated in Fig. 3(a) For comparison, the results of the dynamical systems method 
from section 13.41 are indicated with + symbols. 



Fig. H] also shows the dispersion relation for {E^^D) = (—1,0.1) determined with the 
dynamical systems method in the previous section 13.4) these numerical results are denoted 
with + and can now be compared with the results of the initial value problem from the 
present section. Around the maximum of the curve, the agreement between the numerical 
results of the two very distinct methods is convincing. For larger values of k , the differences 
increase, but the error bars of the initial value problem results increase as well. Furthermore, 
the plotted error bars are an underestimation as they only account for the errors discussed 
above that emerge from the choice of the initial condition and from the time interval of 
evaluation and therefore from possible initial transients and from a possible transition to 
nonlinear behavior. Additional errors can be due to the numerical discretization and time 
stepping of the s themselves. We therefore conclude that the two results agree within the 
numerical error range of the initial value simulations over the whole curve. 
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5 Analytical derivation of asymptotic limits for k <C 1 and 

fc» 1 

Having determined the dispersion relation numerically for different values of electric field £"00 
and diffusion constant D in section El and having tested the correctness of the eigenvalue 
calculation against numerical solutions of the initial value problem in section HI we now will 
analytically derive asymptotic expressions for the dispersion relation for small and large values 
of the wave modes k . It will be shown that these asymptotic limits are 



s(k) 



-Dk 2 , k > 1 



In doing so, we mathematically formalize and generalize the derivation of the small k asymp- 
totic limit that was presented in [3] for the singular limit D = 0, and we correct the result 
proposed in [5]; and we also rigorously derive the large k asymptotic limit, in agreement with 
the form proposed in [5] . 



5.1 Analysis for the asymptotic limit A; <^ 1 

Translation invariance and electrostatic gauge invariance give two explicitly known bounded 
solutions of the linearized system (|3.3p at k = and s = . These are 

u' (0 = W(0,<(0,Po%,K(0,-Eo(Of and e 5 = (0, 0, 0, 0, 1) T . 

Note that e§ is a solution of the linearized system (|3.3p for k = and s arbitrary. 

From the asymptotics of f|3.3|) for k = = s&t£ t = — 00, we see that the only exponentially 
decaying solution at £ = — 00 is given by Uq(£) . This solution is related to the only positive 
eigenvalue fi^_ (see (|3.5p ). For £ — > +00, the solution Uq(£) — > — E^e^, hence this solution is 
not exponentially decaying for £ = +00 . However, it is easy to obtain an explicit exponentially 
decaying solution at £ = +00, this is Ug(£) + E^e^ . 

From the eigenvalues in (|3.5p it follows that for < k <C 1 and < s <C 1 , the three- 
dimensional unstable manifold at £ — > —00 involves one eigenfunction with a fast exponential 
decay (related to the eigenvalue ) and two eigenfunctions with a slow exponential decay 
(related to the eigenvalues k and -4r ) • Similarly, from the eigenvalues in (|3.9p , it follows 
that the two-dimensional stable manifold at £ —* +00 involves one eigenfunction with a fast 
exponential decay (related to the eigenvalue —A* ) and one eigenfunction with a slow 

exponential decay (related to the eigenvalue —A;). Recall that the stable manifold is defined 
as a subset of the full stable manifold to account for the spatial decay condition (|2.1ip . 

We focus on approximating an exponentially decaying solution on the stable manifold. As 
we have seen above, in lowest order, this solution is 

w s (£; £oo, 0, 0) = u' (0 + Eooes + Oik + a). 

To determine the next order, we will use the slow behavior of the asymptotic system and 
write 

w s (£; E^, k, s) = u' (0 + E^ (0, 0, 0, k, 1) e~ fc « + kU^) + sUf^) + + kf). 
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The second term on the right-hand side (-Boo (0, 0, 0, k, 1) e~ k ^ ) is an approximation of the 
slow behavior on the stable manifold, while the other three terms are related to the fast decay. 
Because of the slow decay, the expansion is only valid on a ^-interval with k^ = o(l), hence 
£ shouldn't be too large. 

Substitution of these expressions into the linearized system (|3.3p gives that Uf s and Uf k 
have to satisfy 



p r M(p oo ,0,0)) Uf ta = (%,0, 0,0 

— M(£; Boo, 0, 0)) Uf k = -£?«, (-^^e- fc 5,0,^e-^,0,l- e" fc « 



By analyzing the unperturbed system, we can find a particular solution of the first equation. 
The front solution uq = (tq, ctq, po, Eq, 4>q) satisfies 



J _ ^0 /o 

Po - — 
E 'o = -(°"o - Po) 
= -Bo 

Differentiating this system with respect to B^ gives 

(D£ _ MK;£oo ,„, 0)) ^^(_g, ,^, , )._^ (4 0,^,0,0) 

Hence 

77S / cfo* \~ <9uo . 

tVi „ = — — — — 7— — + a homogeneous solution. 

\dEooJ dEoo 



The asymptotic behavior of J^p- is 

duo 

dEoo ' 



(0,0,0,1,-0. £^00. (5.1) 



The polynomial growth in the 4>k -component will need to be canceled by the behavior of the 
other terms which involve k and hence will give a relation between s and k . 

In fact, the E k -components of the linearized system (I3.3P for any s or k can be expressed 
by the integral equation (|3.4p . For all solutions on the stable manifold, the limit Ek(oo) = 
limf_ >00 Bfc(£) is well-defined, so we can write the E k -components on the stable manifold as 

1 f°° r i 
Kit) = <*e~* - i J [e k ^ + e' k ^\ [p k ( V ) - a k ( V )] dr,, 

In the integral, a k must satisfy the decay condition (|2.11|) on the stable manifold and hence 
will have fast exponential decay. Furthermore, from (|3.3p . it can be seen that pk(£,) = 

c^e^ + -r / e [o'o( ? ?)/o( r ?)Bj! c (r/) — fo(f])ak(f])] drj. As the term inside the integral has 

ho 

fast exponential decay, we get that on the stable manifold C3 = and p k has fast exponential 
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decay too. Thus the integral in the expression for El has fast exponential decay for £ large. 
Since 4> k = ~E' k , we get on the stable manifold for k small and ^-values not too large, say 

1 3 

£ k~i (hence k£ ~ k± ) 

= (k, 1) ( 1 - k£ + O(kVk)) = (k, 1 - kO + O(kVk). 



The exponentially decaying solution on the stable manifold is given by 
w s (Z;E oo ,k,s) = u' (0-s (Jj£-\ 1 ^(0+E oo (0,0,0,k,l)e- k ^ + kUl k (0+O((s+k) 2 ), 

and the arguments above show that the order k contribution in the (E k , 4>k) -components 
is given fully by i?oo(0, 0, 0, k, 1) e~ k ^ and that kU[ k (£) does not contribute to those com- 
ponents at this order. So it follows that the polynomial growth in the (j) k -component of 
-§§^-(0 as given by (|5.ip has to be canceled by the ^-component in i?oo(0> 0, 0, k, 1) e~ fe £ , 



dv* 

s = c*k + 0{k 2 ), with c* = E 00 -j— , (5.2) 

and given in Eq. (|2.8p . Equation (|5.2p establishes the small A: limit of the dispersion 
relation s(k) . 

5.2 A physical argument for the k ^ 1 asymptotic limit 



There is also a physical argument for the asymptotic limit (|5.2p that generalizes the calculation 
in section IV. C of reference [3] to nonvanishing D > 0. For k <C 1, the wave length of the 
transversal perturbation 2ir/k is the largest length scale of the problem. It is much larger 
than the inner longitudinal structure of the ionization front. On the length scale 2ir/k, the 
front can therefore be approximated by a moving boundary between ionized and non-ionized 
region at the position 

z f (x,t) = z + v*{E oo )t + 5e ikx+st , (5.3) 

and the local velocity of this perturbed front is 

v{x,t) = d t z f {x,t) = v*{E OQ ) + s5e ikx+st . (5.4) 

The electric field in the non-ionized region is determined by E = — V0, where <fi is the 
solution of the Laplace equation V 2 = together with the boundary conditions; these are 
E — ► E^z for z — > oo fixing the field far ahead of the front and 0(...) = 0(k) ~ making 
the ionization front almost equipotential. (Due to gauge invariance the constant potential 
can be set to zero.) The solution of this problem is 

</>(x,z,t) = -E oo (z-z -v*t)+E oo e' k(z ' z f ) 5e ikx+st + 0{5 2 ), for z>z f , 
E + (x,t) = E^ + kE^ 6 e ikx+st + 0(5 2 ), for z = z f , (5.5) 
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here E + (x, t) = lim^o E(x, Zf + e, t) is the electric field extrapolated onto the boundary from 
the non-ionized side. 

As the perturbation is linear, the front is almost planar 5 <C 2ir/k. Therefore it will 
propagate with the velocity v*(E + ) = \E + \+ 2y/Df(\E+\) of the pi anar front in the 

local field E + . Inserting E + from (|5.5p and expanding about E^ , we get 



v(x,t) =v*(E + (x,t)) =v*(E OQ )+d E v* k Eoo S e ikx+st + 0(6 



(5.6) 



Comparison of (|5.6p and (|5.4p immediately gives the dispersion relation s = c*k + 0(k 2 ) (|5.2[) . 
that generalizes the result s(k) = \E 00 \k + 0{k 2 ) that was derived in [3] for the singular limit 
D = 0. 



5.3 Analysis for the asymptotic limit k ^> 1 

The asymptotic limit for k 3> 1 is derived by a contradiction argument. We will suppose 
that k is large and that s + Dk 2 is positive, but not small, i.e., 



k > 1, s + -DA; 2 > 0, and s + Dk 2 ^ o(l) 



(5.7) 



and show that this does not allow for bounded solutions. With the assumptions above on s 
and k , the dominant contributions in the matrix M on the whole axis £ are 



(o 


s+Dk 2 








o \ 


D 


1 




















s 

V* 




















-k 2 


Vo 








-1 


o ) 



+ 0(1). 



(5.8) 



Here the three entries —1, 1 and are necessary for a nonvanishing determinant. 

We want to use the Roughness Theorem |13| for exponential dichotomies to show that 
for k large and s not close to —Dk 2 , the exponential dichotomy of the constant coefficient 
ODE is close to the exponential dichotomy of the full system. So first we recall the definition 
of an exponential dichotomy, which gives projections on stable or unstable manifolds. 



Definition 1 ( |13| ) Let A be a matrix in W 
&(y) be a solution matrix of the linear system 

du 



■ , u e M. n , and J 



or 



— = A(y)u, y G J. 
dy 



Let 



(5.9) 



The linear system (15. 9h is said to possess an exponential dichotomy on the interval J if there 
exist a projection P and constants K and k s < < k u with the following properties: 

|*(y)P*- 1 ( 2/Q )| < Ke K> ("»), for y > y , y,y G J 
|*(y)(I - P)*- 1 ^)! < Ke« v <*-v°), for y >y,y,y £j 



An extension for PDEs of this definition can be found in [41] . 

The Roughness Theorem for exponential dichotomies states the following. 
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Theorem 2 (Roughness Theorem [13J ) Consider the system 

[Ao + Ai(y)]«, 



du 
dy 



(5.10) 



with Aq G M riXn a hyperbolic matrix and u £ W 1 . Then for all 5q > there exists a 
5\ > such that for all matrix functions Ai : M — > M 7ixra uiii/i || AiH^oo^+^nxn) < <5i , i/ie 
system (|5,10p /ias an exponential dichotomy on R + (and its dichotomy exponents 

and projections So-close to those of ^jj = Aqu (in the L°°(M. + , M nxn ) norm). 

A constant coefficient linear system does not have bounded solutions. So if the exponential 
dichotomy of the linearized system (|3.3p is close to the one of the constant coefficient system 
with the matrix as in (|5.8j) . the linearized system f)3.3f) does not have bounded solutions 
either. We will show that this is the case if s and k satisfy the assumptions (15. 7ft . 

First we introduce some scaling and coordinate transformations. Define the small pa- 
rameter e = r and the scaled spatial variable, the transformation matrix and transformed 
vector 



£;£ = -, T(e) = diag(e, 1, e, e, 1) and w(n) = T(e)w(en). 



Now (13.31) can be written as 



with 



and 



Mi fa; 





= M (e, 


s) +eM 1 (T);E <x ,e) w, 






/o 


1 + — 


\ 






1 








M (e 


») = 








^00 

V 












0-1 






\o 





-1 ) 


D 




D 


-fo £a 

D 


D 

















— e 


/o 

V* 





^ofo 

V* 







e 


1 


















(5.11) 



where £ = en. 



The eigenvalues and eigenvectors of the constant coefficient matrix Mo are 

±1, with eigenvectors w±i = (0,0,0,^1,1); 

±\J 1 + ^f, with eigenvectors w± 2 = ^±-^/l + 1,0,0,0^ ; 

p-, with eigenvector W3 = (0, 0, 1, 0, 0) . 

If \s\ <C \i then the matrix Mo(s, s) is not hyperbolic at e = 0. However, this problem is 
not fundamental as it is known that there is a hyperbolic splitting in the full problem (see 
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section I3.2f) and there is a spectral gap between the positive and negative eigenvalues, even if 
£ is close to zero. The spectral gap disappears if s ~ —Dk 2 (or e 2 s ~ —D) and in this case 
the following arguments will not work. The spectral gap allows us to define a weight function 
which moves the spectrum away from the imaginary axis. To be specific, define 



w(?7) = e uri 'w(ri), with v = < 
Then Vf(rj) satisfies the ODE 



| sgn(s), if |es| < v* 
0, if \es\ > v* 



Mo(e,a) +ul + eMi (??;£oo, e) ) w (5.12) 



and the spectrum of Mo(e, s) + v\ is bounded away from zero for all e small (as long as 
s + Dk 2 is not small). The system 



M (e,s) + vl w (5.13) 



has an exponential dichotomy with projection Pq such that the range of Po is the span of all 
eigenvectors of the negative eigenvalues and the kernel of Po is the span of all eigenvectors 
of the positive eigenvalues. Then Po is the projection on the stable subspace of the linear 
system (j5. 13|) and I — Po is the projection on the unstable subspace. 

Clearly the matrix Mi (77; , e) is uniformly bounded for all 77 and e small. Thus 
applying the Roughness Theorem [2] gives that for e small there is an exponential dichotomy 
for the system (|5. 12|) on R + with projection P| (77) onto the stable subspace such that P| (7?) 
is e -close to Po for all 77 > 0. And similarly, there is an exponential dichotomy for the 
system (|5.12p on M~ with projection P" (77) onto the unstable subspace such that P"(?7) is 
e-close to I — Po for all 77 < 0. Thus the range of P|(0) is e-close to the range of Po and 
the range of P" (0) is e-close to the range of I — Po , so the range of P*(0) and the range of 
P"(0) have only a trivial intersection for e small. 

As the weight function e vr] has been chosen such that no eigenvalue crosses the imaginary 
axis, it affects only the value of the dichotomy exponentials, not their sign, nor the stable 
and unstable manifolds. Hence the stable and unstable manifolds of (15. lip have a trivial 
intersection only. And the same holds for the stable and unstable manifolds of (|3,3p . as the 
only difference between the systems (|5. 11 [) and (|3.3p is a scaling. So it can be concluded that 
the linear system (|3.3p does not have any bounded solutions for e small ( k large) and s not 
close to —Dk 2 . 

If s is close to —Dk 2 , i.e., s = — + o(l)), then the matrix Aq(b, s) has a positive 
and negative eigenvalue of order o(l) and the spectral gap will disappear in the limit e — > 0. 
So the roughness theorem can not be applied anymore and no conclusion about bounded 
solutions can be drawn. 

The arguments above show that only if s = —Dk 2 (l + o(l)), there is a possibility for 
bounded solutions to exist. As the dispersion curve indicates a bounded solution of the 
ODE (13. 3p . this implies that for e near zero, hence k large, the dispersion curve is 

s(k) = -Dk 2 (1 + o(l)) , k^oo. (5.14) 



So far we have used that s + Dk > 0. If one considers the linear system (|3.3p . an edge 
of the continuous spectrum is given by the curve s = -Dk 2 + /(l^ooD- However, one should 
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include the decay condition (|2.1ip . i.e., the scaling (13, 7p . For any < /3 < A* the edge 
becomes s = —Dk 2 + j3 . By taking the limit for f3 — > 0, we see that the curve s = —Dk 2 
is an edge of the continuous spectrum. Thus with the decay condition, either the dispersion 



6 Physically guided fits to the numerical dispersion relations 

In section [3] we have derived dispersion relations for a number of fields and diffusion 
constants D by numerically solving an eigenvalue problem, we have confirmed these calcula- 
tions by a numerical solutions of the initial value problem in section [H and we have derived 
analytical asymptotic limits to these dispersion relations in section [5j This sets the stage for 
comparing the numerical results to the analytical asymptotic limits and for deriving physically 
guided empirical fits to the numerical dispersion curves where the analytical asymptotic limits 
are not applicable. The small /c-data derived in section [3741 and the analysis of section [57T1 
are shown to be consistent in section 16.11 After showing in section 16.21 that a simple cross- 
over formula joining the asymptotic behavior of the small wave numbers with the asymptotic 
behavior of the large wave numbers is not satisfactory, we give a data collapse, empirical fits 
and arguments on relevant scales in section 16.31 

6.1 Testing the small k asymptotic limits 

First the asymptotic relation (|5.2p for small k is tested on the numerical results. Beyond 
the results visible in the plots of figure [2j the numerical dispersion relation for = — 1 
and D = 0.1 was evaluated carefully for small values of k, and the result is shown in a 
double logarithmic plot in Figure [5] that zooms in on the small k behavior. Also plotted is 
the analytical asymptotic limit (|5,2p . The comparison between numerical data and analytical 
asymptotic limit is convincing in the range of small k . 

For the other values of E^ and D presented in figured the dispersion relation s(k) again 
is fitted very well by the asymptotic limit (|5.2p for small values of k . This will be illustrated 
in Figure [6] below. 

6.2 Testing both asymptotic limits 

It is quite suggestive to join the small k asymptotic limit (j5.2|) with the large k asymptotic 
limit (|5.14p into one cross-over formula 



A formula similar to s(k) = c*k — Dk 2 was suggested in JS], but with a different prefactor 
instead of c* . However, we now will confirm once more the correctness of the prefactor c* , 
and we will show that the large k asymptotic limit is not yet applicable in the range of 
positive growth rates s(k) . 

If (|6.ip holds, then the dispersion relation for the rescaled variable S = Ds/c* 2 as a 
function of the rescaled wave number k = Dk/c* becomes S = k — k 2 . Therefore the formula 
(16. ip can easily be tested on the numerical data from figure [2] by plotting them in rescaled 
variables S and k with appropriate values for D and c*(-£7oo, D) for each curve. The result is 



curve satisfies s(k) > —Dk 2 



or it ends at the continuous spectrum. 




(6.1) 
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Figure 5: A log-log plot of the dispersion curve s(k) for E = — 1 and D = 0.1 illustrates 
the behavior for small k. Red squares and solid line: Data from the numerical evaluation of 
the eigenvalue problem. Blue dashed line: Analytical asymptotic limit log s = log c* + log k 
according to (|5.2p . 

shown in Fig. [6l together with the parabola k — k 2 . The plot illustrates that the asymptotic 
limit (|5.2p indeed is a very good fit to all data for small k . 

For larger k, the curves differ quantitatively. In particular, S vanishes for k between 
0.014 and 0.035 for the numerical dispersion curves while the formula (|6.ip predicts this to 
happen for n = 1 . Also the maximum of the dispersion curve cS max is never higher than 0.0027 
for the numerical data while formula (|6.ip predicts 0.25. Of course, this is not in contradiction 
with the analytical results in Section 15.31 for large k . Rather it says that the positive part of 
the dispersion curve lies completely in the range of small k, where the asymptotic limit for 
large k is not applicable. We conclude that cross-over formula (16. ip is not an appropriate fit 
for the numerically derived dispersion relations. 

6.3 Data collapse, relevant length scales, empirical fits and conjectures 

We finish this section with a data collapse and arguments on relevant scales that guide em- 
pirical fits. 

6.3.1 Data collapse 

First, we investigate whether the numerical data for different E^ and D can be collapsed 
onto one curve. This is done by determining the maximum of the dispersion curve s max and 
the wave number ko where the growth rate vanishes, s(ko) = 0, from the numerical data for 
each pair (Eo^D) . In figure[7]all curves are plotted as s/s max and k/ko with their respective 
s max and ko . The plot shows that the curve shapes are very similar, but they do not coincide 
completely. For example, there seems to be a small drift in the position of the maximum. 
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Figure 6: Labeled curves: The numerical dispersion curves from Figure [5] plotted as S = 
Ds/c* 2 over k = Dk/c*. Dotted line on the left: the parabola k — k 2 that would be 
predicted by (|6.ip as far as it fits into the plotted region. 



6.3.2 Relevant length scales and the D = case 

In a second step, we investigate which physical or mathematical mechanisms can suppress 
the growth rate s(k) for much smaller values of k than suggested by the large k asymptotic 
limit (|5.14p . In a first overview, there are three length scales in the problem. The transversal 
perturbation is characterized by its wave length 2n/k. In the longitudinal direction, the front 
is characterized by two length scales, the electric screening length t a and the diffusion length 
i D , cf. Fig.ffl 



l a = v and £ D = — = J — ^ (6.2) 

For vanishing diffusion D, the diffusion length vanishes, and the screening length l a has to 
be compared to the wave length of the perturbation; in [3] it was shown that it determines 
the cross-over from the small to the large k asymptotic limit of the dispersion relation: 

s(k) = T\u <<k r ^D = 0, where k a = -L ( 6 .3) 

V 1 I |£oo| k a , for \k\ >/c Q ; 2£ a v ' 



The actual curve for D = and = — 1 is given in figure 2(b) , where we remark that the 
form 

8(k) = ^ m 

for positive real p reproduces the asymptotic limits (|6.3p . but does not fit the full numerical 
curve for = — 1 satisfactorily for any power p. The functional form of (16. 4h will serve 
below as an inspiration for our empirical fits for D > 0. 
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Figure 7: The numerical dispersion curves from Figure [2] plotted as s/s max over k/kq] here 
s max = maxfc s(k) and ko > with s(ko) = are determined from the respective curve. 



Searching for why (|6.4p does not properly fit the data, one realizes that there are actually 
two different definitions of the screening length possible: 

j f =*(E 0O ) and ^ = X~= f ]EO ° ] ^^ + 0(VD), (6.5) 

with A - from (I2.17p . The dimensions of both quantities are the same, and they approach each 
other if a is constant in a large part of the integration interval [0, |-Eoo|] j this is the case with 
the Townsend approximation a(x) = e~ x l x for l-E^I ^ 1- Otherwise, £+ characterizes the 
slopes of the fields near the discontinuity of a [3], while £~ characterizes the decay (I2.14p of 
the fields far behind the front for £ — * — oo . The analysis in [3] shows that linear perturbations 
with wave numbers k 3> 1 couple to the inner local structure of the front and are dominated 
by £+ , while smaller k could couple to the larger spatial structure characterized by £~ , this 
conjecture will be tested on the numerical data below and asks for future analysis. 



6.3.3 Scales and fits for D > 

When diffusion is included, the diffusion length £o emerges as another length scale in the 
front. As illustrated in figure [H instead of the discontinuous electron density in the front for 
D = 0, a diffusive layer of width £rj = 1/A* (I2.13|) builds up in the leading edge. While 



D increases, the dispersion relation decreases as shown in figure 2(b) As the diffusive layer 
is the main new feature of the front for D > 0, it is plausible that the different behavior 
of s is created within this boundary layer. The physical mechanism is that diffusion can 
smear perturbations of short wave length out, hence suppressing their growth. This process 
mainly takes place in the diffusive layer because gradients are largest in this region. This 
idea has inspired an attempt in [5] to calculate s(k) by local analysis within the diffusion 
layer. In principle, such an approach combined with proper matched asymptotic expansions 
could work. However, the calculation in [5] was intrinsically inconsistent [18j . disagrees with 
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our asymptotic limit for small k and therefore fits the numerical results even worse than 
formula ||0]L cf. Fig. [3 

We have tested whether the diffusion length £d = 1/A* plays a role in the dispersion 
relation by plotting the numerical data from Fig. [2] this time for the rescaled variables s = 
s/(c*A*) over 8. = k/A*. The result is shown in Figured) It shows that the numerical 
dispersion curves are well approximated by 



s(k) m c*k + 0(k 2 ), for/fc^O, 
k sa A*/4, where s(k ) = 0, 



(6.6) 
(6.7) 



The numerical evidence from Fig. [8] summarized in (|6.7jl together with the physical explana- 
tion above suggest the following conjecture. 

Conjecture 1 T/ie largest unstable wave number of the Laplacian instability is proportional 
to the inverse diffusion length. 

We remark that the data gives hy ~ 1/(4£d) , while the cross-over formula (|6.ip would suggest 
that ko « £ a /£ 2 D , highlighting again its inadequacy for intermediate k values. Figure [8] also 
shows that the value of the wave number for which the maximum of the dispersion curve is 
attained, lies in the range of k max = 0.22 ko to 0.30 ko . 
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(b) Using a = 3c*/f(\E oa \) in formula (gTS 



(a) Using a = 3/a(\Eao\) in formula (|6.8 



Figure 8: The data of Figure [2] plotted as s = s/(A*c*) over R = k/A* . The lines are given 
by the empirical formula (16.80 . 



The data in Figure [8] suggest an empirical formula of the form (for s > ) 
c*k / 4k\ 8. 



s(k) 



l + ak 



1 



4fc 
A* 



or 



1 + aA*R 



(6.8) 



where the parameter a will depend on the external parameters D and E^. The factor c*k 
creates the correct asymptotic limit (j6.6j) for fe < 1. The factor (1 — 4£;/A*) creates the 
non-trivial zero of the dispersion relation at ko (|6.7p . The form of the numerator is inspired 
by (|6.4[) . and the proper asymptotic limit (16.3H for large fc and D = would be reached for 
a = 2£+ + C(v / ^) • Obviously, the empirical formula (|6.8|) is not valid in the asymptotic range 
k 1 where s < and where the asymptotic behavior is given by s ~ — Z?/c 2 . 
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The functional form of formula (|6.8H is supported by the following observation. If one cal- 
culates the maximum (^ m ax>s max ) of (|6.8p . it follows that ^ max /s m ax = 1/4, independently 
of the value of a. (The number 1/4 directly stems from the factor 4 in (1 — 4^).) This 
relation indeed fits the numerical curves quite well, therefore the factor 4 is supported twice 
independently. Relevant numerical data for this and other fits is collected in Table [TJ 

The value for a is less obvious. The empirical formula (|6,8p gives the following relation 
between aA* and the maximum of curve 

1 — 8i£ ma x ^ ~~ 4\/ ^max 

4-^max Smax 

The empirical values for those quotients are given in Table [TJ 

The limit of D = and I; > 1 suggests a = 2£+ + 0(\^D) , but the fit is unconvincing 
(cf. also the discussion for D = above). However, we found that a = 3£+ fits the data 
reasonably well. Formula (|6,8p with this value of a together with the numerical data are 
visualized in Figure [8]^a). The fit is quite good for the lower two curves, the upper two 
display some discrepancies. The main problem is the value of s max with a relative error 
between 2% and 24%, while the position of ^ m ax has a relative error as low as 0.5% to 7%. 
As /(|£oo|)/c* = a(|-Eoo|) + 0{Jl5), another possible fit is a = 3c* / /(|£oo|) , it is displayed 
in Figure E^b). The fit is quite good for the upper two curves, but now the fit has some 
discrepancies for the lower two curves. And the value of ^ max has a larger error between 2% 
and 10% while the position of s max has a much smaller error of only 1% to 10%. Obviously, 
these observations ask for further analytical investigation. Note finally the striking relation 
between A~ = and the value of /c max for larger values of the electric field in Table [TJ 

As a basis for future work, all characteristic numerical data is collected in this table. 

7 Conclusion and outlook 

In this paper, we have found dispersion curves for negative streamer ionization fronts by 
numerically solving an eigenvalue problem, we have verified this prediction on the numerical 
solution of an initial value problem, we have derived analytical expressions for the asymptotics 
of the curve for large and small wave lengths, and we have presented a physically motivated 
fit formula to the numerical curves for intermediate wave lengths. The investigation is of 
interest for two reasons: because pulled fronts like these ones are mathematically challenging 
to investigate, and because explicit predictions on the linear stability of ionization fronts help 
to interpret numerical and experimental observations of propagating and branching streamer 
discharges. 

The ionization front is a pulled front, i.e., the front is part of a family of traveling waves, 
which propagate into a temporally unstable steady state. For the dynamics with one spatial 
variable, most traveling waves in this family are attractors only for waves with exactly the 
same asymptotic decay profile. The exception is the pulled front, which has the steepest 
decay of all waves in the family and is an attractor for waves with a sufficiently fast decay 
(therefore excluding the slower decay rates for the other traveling waves in the family). The 
instability of the state ahead of the front and the related spatial decay condition imply that 
only a submanifold of the stable manifold in the transverse instability problem is relevant for 
the transverse instability analysis. This submanifold is identified by introducing a weighted 
solution space that excludes solutions with a too slow decay rate. We have integrated the 
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Table 1: Upper block: characteristic numerical data of the dispersion relations in Figure [2] 
with errors. Middle block: characteristic scales of the planar front according to analysis. 
Lower block: relevant ratios of numerical and analytical scales as used for the derivation of 
the empirical formula. 



{Eoo, D) 


(-1, 0.01) 


(-1, 0.1) 


(-2, 0.1) 


(-5, 0.1) 


(-10, 0.1) 


^max 


0.080(1) 


0.05190(2) 


0.1695(15) 


0.647(1) 


1.6305(15) 


h 

"inax 


0.35(4) 


0.144(1) 


0.25(4) 


0.45(4) 


0.60(4) 


k 


1.575(15) 


0.4825(5) 


0.875(25) 


1.595(1) 


2.397(1) 




4.56(56) 


3.35(3) 


3.60(65) 


3.57(33) 


4.01(27) 


V* 


1.12 


1.38 


2.70 


6.28 


11.9 


c* 


1.12 


1.38 


2.52 


5.77 


11.0 


a* = v7(|£ool)/£> 


6.07 


1.92 


3.48 


6.40 


9.51 


a(l^ool) = i/e+ 


0.37 


0.37 


0.61 


0.82 


0.90 




0.148 


0.144 


0.638 


2.832 


7.169 


A- = 1/ia 


0.13 


0.10 


0.23 


0.45 


0.60 


SAVadEool) 


49.5 


15.6 


17.2 


23.4 


31.5 




55.5 


21.6 


21.7 


27.0 


34.8 


&o = k /A* 


0.260(3) 


0.252(1) 


0.251(5) 


0.249(1) 


0.252(1) 


Smax — Smax/ ip A ) 


0.0118(2) 


0.0196(1) 


0.0193(2) 


0.0175(1) 


0.0155(1) 




1.9(2) 


0.78(1) 


0.8(2) 


1.1(1) 


1.3(1) 


■^maxAmax 


0.28(7) 


0.288(4) 


0.27(9) 


0.28(5) 


0.26(4) 


(l-8^ max )/«J 


40(6) 


17.7(1) 


21(2) 


22(2) 


31(2) 


(1 — 4-y/s max )/s max 


48.1(4) 


22.5(1) 


23.0(1) 


26.8(1) 


32.3(1) 



relevant stable submanifold and unstable manifold numerically with a dynamical systems 
method to calculate the dispersion curve. This method of finding the dispersion curve does 
not use any details of the streamer model, except that it has a pulled front. The definition of 
a submanifold of the stable manifold and the subsequent numerical integration of this stable 
submanifold and the unstable manifold are ideas that can be applied to pulled fronts in other 
systems, too. 

It is interesting to see that the band of unstable wave numbers seems to be limited by a 
multiple of the decay rate A* that characterizes the leading edge of the pulled ionization front; 
though the evidence up to now is only numerical. As such behavior is physically reasonable, 
the next step would be to derive it analytically, e.g., by a local analysis in the diffusive layer 
and matched asymptotic expansions. Such an expansion could be based on the limiting case 
where the diffusion length Id = 1/A* is much smaller than the screening length £ a . 

The calculated dispersion curves also contribute to understanding the stability of actual 
streamers. Two- and three-dimensional time dependent simulations [U I4*2l [36l [37} [33] of the 
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streamer model introduced in section 2 show them to become unstable and branch. Can the 
unstable wave lengths of this branching be related to the unstable band of wave lengths of the 
present calculation? Furthermore, if the inner front structure is approximated by a moving 
boundary |35|. I19j . how is the calculated dispersion relation of transversal perturbations to be 
taken into account? It will also be interesting to see whether the dispersion relation calculated 
for the present fluid model is also applicable to the corresponding particle model [32J. 

Finally, we mention that the extension of the streamer model with photo-ionization as an 
additional reaction term [33] in composed gases like air requires an extension of the present 
analysis as nonlocal interaction terms play a role. 

A Matrices in exterior algebra spaces 

In this appendix, we give explicit expressions for the matrices acting on the exterior 

algebra space /\'(C 5 ) for I = 2,3. Let ex,..., 65 be the standard basis for C 5 . Then an 
induced basis on /\ 2 (C 5 ) is given by 

ai = ei A e 2 , a 2 = ei A e 3 , a 3 = ei A e 4 , a 4 = e x A e 5 , a 5 = e 2 A e 3 , 
a 6 = e 2 A e 4 , a 7 = e 2 A e 5 , a 8 = e 3 A e 4 , a 9 = e 3 A e 5 , a 10 = e 4 A e 5 . 

The matrix M^ 2 ) : A 2 (C 5 ) -» A 2 (C 5 ) can be associated with a complex 10 x 10 matrix 
with entries such that 

10 

M^a^^M^a,-, i, j = 1, . . . , 10 , (A.l) 

i=i 

where, for any decomposable x = xi A x 2 G /\ 2 (C 5 ) , M^x := Mxi A x 2 + xi A Mx 2 . Let 
M be an arbitrary 5x5 matrix with complex entries, 





mn 


m 13 


mu 




m 21 


m 2 2 


"1-23 


m 24 


m 2 5 


m 31 


m.32 


"133 


m 34 


m 3 5 


m 4i 


m 42 


m 43 


m 44 


m 45 


\m 5 i 


m 5 2 


"153 


m 54 


^55/ 



(A.2) 



then M^ 2 ) takes the explicit form 



M( 2 ) = 





"123 


m 24 


m.25 


-77113 


-mu 


-77115 








" 


m 32 


d\3 


m 34 


m 3 5 


"112 








-77114 


-77115 





m 42 


m 43 


du 


m 45 





"112 





mi3 





-77115 


"*52 


"*53 


"2-54 


diB 








mi 2 





mi 3 


"114 


-m 3 i 


m 2 i 








d23 


m 34 




-m 24 


-m 2 5 





-m 44 





m 2 i 





771 43 


^24 


77145 


"123 





-"I25 


-JTO51 








m 2 i 


"I53 


m 54 


^25 





"123 


m 24 





-77141 


"131 





-m 42 


"i 32 





CZ34 


m 45 


-"135 





-77t 5 l 





m 31 


-m 52 





m.32 


77154 


45 


m 34 








-77151 


77141 





-"152 


m 42 


-m 53 


m 43 


^45 . 



where du = ma+mjj. 
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In a similar way, the matrix : /\ 3 (C 5 ) -> A 3 (C 5 ) can be associated with a complex 

10 x 10 matrix. First we define an induced basis on /\ 3 (C 5 ) by 

bi = ei A e 2 A e 3 , b 2 = ei A e 2 A e 4 , b 3 = ei A e 2 A e 5 , b 4 = e 1 A e 3 A e 4 , 
b 5 = ei A e 3 A e 5 , b 6 = ei A e 4 A e 5 , b 7 = e 2 A e 3 A e 4 , b 8 = e 2 A e 3 A e 5 , 
b 9 = e 2 A e 4 A e 5 , bi = e 3 A e 4 A e 5 . 

The matrix for M^ 3 -* £ C 10xl ° has entries such that 

10 

M^b^^Mfbj, i,j = l,...,W, (A.3) 
i=i 

where, for any decomposable x = x 4 A x 2 A x 3 £ /\ 3 (C 5 ) , M^x := Mx 4 A x 2 A x 3 + x 4 A 
Mx 2 A x 3 + xi A x 2 A Mx 3 . If M is given by (IA.2j) . then M^ 3 ) takes the explicit form 



m( 3 ) = 



^123 


m 43 


"153 


-"242 


-"252 





"141 


77251 





" 


"2-34 


d\2A 


"254 


"232 





-"252 


-"131 





"251 





m 35 


m 45 


^125 





"232 


"242 





-"131 


-"141 





-m 24 


"123 





^134 


"254 


-"153 


"121 








m 5 i 


-m-2,5 





"223 


"245 


^135 


"143 





"121 





-"141 





-?"25 


"224 


-"135 


"234 


^145 








"221 


"131 


mu 


-"2-13 





"212 








<^234 


"154 


-"153 


"252 


mis 





-"213 





mi2 





"145 


^235 


"243 


-"142 





"2-15 


-"214 








"112 


-"135 


"134 


^245 


"232 











"215 


-"214 


"113 


"225 


-"124 


"123 


^345 . 



where d ijk = mu+mjj + m kk . 
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